#####
#####
###     Creating Figure 3
#####
#####

library(foreign)
library(lattice)

###
###   Read value coordinates from file, 
###   "value coords, with names.txt"
###

val.coords <- read.table(file.choose(), header = T)

###
###   Religious denomination. Read data
###   from Stata dataset, "relig.dta"
###

relig <- read.dta(file.choose())

###   Create new variable with
###   character values for religious
###   denominations
###

relig$rel[relig$religion == 1] <- "Prot"

relig$rel[relig$religion == 2] <- "Cath"

relig$rel[relig$religion == 3] <- "Jew"

relig$rel[relig$religion == 4] <- "Muslim"

relig$rel[relig$religion == 5] <- "None"

relig$rel[relig$religion == 6] <- "Other"

relig$rel[relig$religion == 7] <- "Christian"

relig$rel[relig$rel == "Muslim"] <- "Other"

###
###   Calculate mean vectors and mean resultant
###   vector lengths for religious denominations
###

mean.rel.d1 <- tapply(relig$dim1, relig$rel, mean)

mean.rel.d2 <- tapply(relig$dim2, relig$rel, mean)

rel.mean.lengths <- (mean.rel.d1^2 + mean.rel.d2^2)^.5

###
###   Plot mean vectors for religious denominations
###   along with value points.
###
###   THIS CREATES FIGURE 3A
###

xyplot(dim2 ~ dim1, data = relig,
       aspect = 1,
       panel = function (x, y) {
         panel.xyplot(val.coords$dim1[-3], val.coords$dim2[-3], pch = 4, 
                    cex = .5, col = "black")
         panel.text(val.coords$dim1[-3], val.coords$dim2[-3], 
                    labels = c("Freedom", "Equality", "Economic\nsecurity",  
                               "Morality", "Individualism", "Social\norder",       
                               "Patriotism")[-3],
                    pos = 4, cex = .5)
         panel.segments(rep(0, 2), rep(0, 2), 
                        mean.rel.d1, mean.rel.d2, lwd = 1.5)
         panel.text(mean.rel.d1[2:5], mean.rel.d2[2:5], 
                    labels = names(mean.rel.d2)[2:5], pos = 4)
         panel.text(mean.rel.d1[1], mean.rel.d2[1], 
                    labels = names(mean.rel.d2)[1], adj = c(.5, 1.1))
         panel.text(mean.rel.d1[6], mean.rel.d2[6], 
                    labels = names(mean.rel.d2)[6], adj = c(.7, 1.1))
       },
       xlab = "",
       ylab = "",
       scales = list(draw = FALSE)
)

###
###   ANAVA on religious denomination
###

meand1 <- mean(relig$dim1)

meand2 <- mean(relig$dim2)

mean.res.length <- (meand1^2 + meand2^2)^.5

res.length <- length(relig$dim1) * mean.res.length

rel.between <- sum(rel.mean.lengths * as.vector(table(relig$rel))) - res.length

rel.within <- length(relig$dim1) - sum(rel.mean.lengths * as.vector(table(relig$rel)))

dfnum <- 6

dfdenom <- length(race$dim1) - 7

rel.ms.between <- rel.between / dfnum

rel.ms.within <- rel.within / dfdenom

F.rel <- rel.ms.between / rel.ms.within

obs.prob.rel <- 1 - pf(F.rel, dfnum, dfdenom)


###
###   Religious commitment, using information
###   from relig data frame
###

###
###   Create a four-category version
###   of the religious commitment variable
###

relig$com[relig$commit >= 0 & relig$commit <= 2] <- "Very low"

relig$com[relig$commit >= 3 & relig$commit <= 5] <- "Mod low"

relig$com[relig$commit >= 6 & relig$commit <= 8] <- "Mod high"

relig$com[relig$commit >= 9 & relig$commit <= 12] <- "Very high"

###
###   Calculate mean vectors for the four
###   categories of religious commitment, along
###   with the mean resultant vector lengths
###   for each commitment level.
###

mean.com.d1 <- tapply(relig$dim1, relig$com, mean)

mean.com.d2 <- tapply(relig$dim2, relig$com, mean)

com.mean.lengths <- (mean.com.d1^2 + mean.com.d2^2)^.5

###
###   Plot mean vectors for four levels of religious
###   commitment, along with value points.
###
###   THIS CREATES FIGURE 3B
###

xyplot(dim2 ~ dim1, data = relig,
       aspect = 1,
       panel = function (x, y) {
         panel.xyplot(val.coords$dim1[-6], val.coords$dim2[-6], pch = 4, 
                    cex = .5, col = "black")
         panel.text(val.coords$dim1[-6], val.coords$dim2[-6], 
                    labels = c("Freedom", "Equality", "Economic\nsecurity",  
                               "Morality", "Individualism", "Social\norder",       
                               "Patriotism")[-6],
                    pos = 4, cex = .5)
         panel.segments(rep(0, 2), rep(0, 2), 
                        mean.com.d1, mean.com.d2, lwd = 1.5)
         panel.text(mean.com.d1[c(2,4)], mean.com.d2[c(2,4)], 
                    labels = names(mean.com.d2)[c(2,4)], pos = 4)
         panel.text(mean.com.d1[c(1,3)], mean.com.d2[c(1,3)], 
                    labels = names(mean.com.d2)[c(1,3)], pos = 1)
       },
       xlab = "",
       ylab = "",
       scales = list(draw = FALSE)
)

